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Abstract. - Networks of caustics can occur in the distribution of particles suspended in a 
randomly moving gas. These can facilitate coagulation of particles by bringing them into close 
proximity, even in cases where the trajectories do not coalesce. The long-time morphology of 
these caustic patterns depends upon the Lyapunov exponents Ai, A2 of the suspended particles, 
as well as the rate J at which particles encounter caustics. We develop a theory determining 
the quantities J, Ai, A2 from the statistical properties of the gas flow, in the limit of short 
correlation times. 



Aerosols are usually unstable systems, in that the suspended particles eventually coagulate. 
Understanding the process giving rise to this coagulation, and determining the time scale over 
which it occurs are important questions in describing any aerosol system. If the gas phase does 
not have macroscopic motion, the coagulation may be effected by diffusion of the suspended 
particles, or (if the suspended particles are of a volatile substance) by Ostwald ripening. The 
coagulation process can be greatly accelerated if the aerosol undergoes macroscopic internal 
motion. Ultrasound, for example, has been used to accelerate coagulation in aerosols [1]. 
Turbulent flow could also play a role in the coagulation of suspended particles; this could be 
relevant in the coalescence of visible moisture into rain droplets [2] . 

If suspended particles are simply advected in an incompressible flow, their density remains 
constant. Inertial effects are therefore required for coagulation, unless the flow exhibits sig- 
nificant compressibility. In earlier work [3,4] we discussed the motion of inertial particles in a 
random velocity field. We showed that there is a phase where the trajectories of the particles 
coalesce, so that arbitrarily small particles coagulate. In the limit where the correlation time 
r of the flow approaches zero, this path-coalescing phase only exists when the velocity field is 
predominantly potential flow (such as the flow due to sound waves) [4]. 

Turbulent fluid flow is, however, expected to be predominantly solenoidal, and it is of inter- 
est to find alternative mechanisms of coagulation which operate outside the path-coalescence 
phase. Here we describe an alternative mechanism facilitating coagulation, illustrated in 
Fig. ^ we show the distribution of particles suspended in a randomly moving gas (the equa- 
tions of motion and statistics of the flow field are given by eqns. ([3) to below). The large 

© EDP Sciences 



2 



EUROPHYSICS LETTERS 




Fig. 1 - Distribution of inertial particles suspended in a randomly moving fluid (blue corresponds to 
lowest density, yellow to highest). The initial distribution is a random scatter. The large panel shows 
caustics at short time. Panels (a)-(c) show the long-time behaviour. In all cases, the region is the unit 
square, the mean particle density is 2.5 x 10^, m — 1, and there is potential flow (with parameters 
5 = 0.03, a = 0.01, St = 0.05, see text). Main panel: 7 = 0.53, t = 5, (a): 7 = 1.18, t = 500, (b): 
7 = 0.72, t = 125, (c): 7 = 0.21, t = 125. The three cases correspond to: (a) A2 < Ai < 0, (b) Ai > 0, 
Ai + A2 < 0, and (c) Ai > 0, Ai + A2 > 0, see text. 



panel shows the distribution of particles after a short time, starting from a random scatter 
with uniform density. The particles cluster onto a network of caustic lines, analogous to the 
networks of optical caustics that can be seen on the bottom of a swimming pool [5] . This phe- 
nomenon is a new mechanism by which aerosol particles are brought into close proximity. The 
remaining parts of Fig.^show the distribution of particles after a long time, in three different 
cases: part (a) shows the path-coalescence phase where the trajectories condense onto points. 
Parts (b) and (c) show two cases where there is no path coalescence, but a steady state with 
significant inhomogeneities of density due to caustics: these have very different morphologies, 
depending on the parameter values, as we shall show. 

Fig. n is surprising because it is be expected that random movement of uniformly dis- 
tributed particles would leave the distribution uniform. The following questions naturally 
arise. First, why do the particle trajectories coalesce into points in Fig. QJa)? This phe- 
nomenon was first noted in [6] and subsequently analysed in detail in [3, 4] (c.f. also the 
theory developed at the end of this paper). Secondly, why and how does the caustic pattern 
develop? We will explain why the caustics form, and relate the pattern to optical caustic 
networks in swimming pools. Third, why do the morphologies of the structures seen in the 
steady state at long times differ? We will argue that the different morphologies of the caustic 
patterns are related to three parameters: the rate J at which any given particle crosses a 
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caustic and the Lyapunov exponents Ai and A2 of the flow (with Ai > A2). We present a 
new theory for these parameters in the Umit where the correlation time of the random flow is 
short. Finally, it is natural to ask how the density fluctuations are quantifled. We will show 
that the probability distribution function for the particle density has two regimes of algebraic 
decay. The formation of caustics causes particles to pass through regions of greatly increased 
density, where coagulation by contact interaction is much more likely to occur. In this paper 
we confine ourselves to describing the caustics, and do not model the coagulation process. 

There is a large literature on light-intensity fluctuations in twinkling starlight, and other 
cases of propagation through random media. Such problems are different because the finite 
wavelength of the light causes the fluctuation distribution due to catastrophes to be cutoff at 
high intensities [5]. Caustics also occur in Hamiltonian dynamics, and their observation in 
electron flows is a topic of current research [7] , where the statistics of caustic formation are of 
interest [8]. Our system differs from the Hamiltonian case: the drag force on the particles is 
assumed to be given by Stokes's law. Because of the Stokes damping, the density correlation 
function reaches a statistically stationary state, as opposed to the Hamiltonian case. Although 
our results are described in the context of aerosols, the same principles apply to suspensions 
of particles in liquids and to tracer particles used to investigate turbulent flows. Because our 
objective is to explain the theoretical principles as clearly as possible, we confine discussion 
to one or two spatial dimensions. 

We consider small spherical non-interacting particles of mass m, radius a, in a random 
flow with velocity fleld u{r, t) and viscosity rj. We assume that the drag force on the particles 
is given by Stokes's law. Neglecting displaced-mass effects, the equation of motion is 

r — — 7(7" — It) (1) 

where 7 = Girrja/m and r is the particle position. This model is widely used in studies 
of particles suspended in turbulent fluids (see [9], [10] for descriptions of recent numerical 
investigations of this model), and it is surprising that the occurrence of caustics was not 
noted earlier. The random driving force on the particles, / — rwyu, is conveniently described 
by two scalar potentials and ip: 

f{r,t) ^\7(j)ir,t) + \7 An3^{r,t) (2) 

(where fi^ is a unit vector perpendicular to the plane): the scalar fields </> and "0 generate, 
respectively, the potential and solenoidal components of the flow. Here we assume that (f> 
and tp are independent, with ((f)) = {ijj) — and {ijj^) = a^(0^) for some constant a (angular 
brackets denote expectation values). Also, we assume that the correlation function of ip has 
the same form as that of 0: 

{(P{r + RM+t)(j,{r,to))^C{R,t), (3) 

where R = and C{R,t) has correlation length ^ and correlation time r. The theory is 
readily extended to more general statistics. 

Our two-dimensional numerical simulations are performed in the limit of small t, using a 
model discretised in time with a small time step St ^ r; the impulse 

l'{n+l)St 

fn{r)= / dt' f{rt>,t') (4) 

JnSt 

at time nSt is taken to be of the form |2l in terms of scalar fields ipnix) and "ipnif) satisfying 
{Mr)(pn'ir'))=a^eSteM-\r-r'\y2e)6nn' and (^„(r)^„, (r')) = a'(<^„(r)^„' (^'))- 
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The origin of the caustics is most easily understood by a one-dimensional example. Fig. [21 
(a) is a schematic plot of the momentum p of a particle at position x. At time t = 0, this is 
a single-valued function, but at a later time this function may develop a pair of folds because 
the faster particles overtake the slower ones. If the density of particles is smoothly distributed 
along the line in phase-space, the projected density in coordinate space, shown in part (b), 
is singular at a pair of caustics, which are the projections of each fold. We see that the 
caustics are created in pairs, and that initially there is a high density of particles between the 
caustic lines, as well as a divergent density on the caustic itself. Fig. |2 (c) shows a segment 
of the actual caustic pattern. The caustic lines cannot abruptly end, but caustics can join 
at cusps [5], some of which are visible in the figure. Some regions of high density do not 
have an associated pair of caustic lines. This is because the phase-space manifold is nearly 
perpendicular to the coordinate-space plane, but it has not yet folded over. At short times 
our caustic networks are equivalent to the lines of bright light observed on the bottom of a 
swimming pool on a sunny day when the water surface is disturbed. The morphology of the 
swimming-pool caustics is discussed by Berry [5] . 

There are important differences between optical realisations of caustic patterns behind 
randomly refracting screens and the caustics in particle density. The former become more 
confused as the distance from the screen is increased, whereas our patterns reach a statistically 
steady state at large times. Although more caustics are created by folding, the damping term 
(with coefficient 7 in eq. ^ ) implies that momentum differences decrease, so that the caustics 
become progressively more tightly folded. 

The examples in Fig. ^b) and (c) show that the morphology of the patterns in this steady 
state depends upon the statistics of the random velocity field. It is desirable to understand 
the differences between these figures. 

We consider the behaviour of three particles, a reference particle and two nearby ones 
separated by Sr and 6r' . The evolution of these small increments is described by the Lyapunov 
exponents, Ai and A2, characterising the time-dependence of the length X{t) — \dr{t)\ of 
one of the separation vectors, and of the area of a parallelogram spanned by two vectors. 




Fig. 2 - (a) Particles are distributed on a phase-space manifold, shown here as a phase curve in a one- 
dimensional section. The phase curve develops folds, (b) The particle density diverges on caustics, the 
projections of the folds. The caustics are created in pairs, with a high density of particles between 
each pair, (c) The particle distribution is shown in red, and the corresponding caustic curves are 
plotted as blue lines. The parameters are the same as for the largest panel of Fig. 
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Fig. 3 - (a): The Lyapunov exponents Ai and A2 determine the particle-density fluctuations, (b): 
Schematic illustration of the model for the folding process leading to eq. (c): Density fluctuations 
in a one- dimensional model — 0.05, m = 1, r = 1, cr = 10~^, 7 = 0.04, averaged over 25 times 
ranging from 40r to lO^r): P{p) shows algebraic decay, with a region of decay followed by a 
tail, (d): Lyapunov exponents Ai and A2 for F = 1/3; simulations of Langevin equations (lines) and 
integration of equation of motion (o). (e): Corresponding rate of crossing caustics. 



A{t) = \6r{t) A Sr'{t)\. These are defined by 



Ai = lim —log 

t^oo t 



X{t) 



X{0) 



Ai + A2= lim -logp 

t^OQ t 



A{t) 



(5) 



Now we use the Lyapunov exponents to gain an understanding of the patterns in Fig. ^ 
When Ai < 0, almost all infinitesimal line segments contract to a point, and the trajectories of 
particles coalesce as shown in Fig.^a). This was the principle used to explain the coalescence 
transition, discussed in [3] and [4]. Next consider what happens if Ai is positive, so that the 
trajectories do not coalesce into points. In this case we consider a small element of area on a 
caustic line, extended in the direction of the caustic [Fig.EJa)]. This element is expected to be 
stretched along the direction of the caustic, because we assume Ai > 0. If A2 is also positive, 
the line width of the element will also expand in the direction perpendicular to the caustic, 
and the density fluctuation resulting from the caustic becomes weaker as time proceeds: this 
is illustrated by Fig.^c). If A2 < 0, the line element contracts in the direction perpendicular 
to the caustic, and if in addition Ai + A2 < 0, the concentration of particles on the caustic 
increases. This case is illustrated by Fig. ^b), which shows very narrow caustic lines. The 
increase in particle density on the caustic lines cannot proceed without limit: the caustic line 
stretches and folds, until the particles aligned along the caustic have become indistinguishable 
from points randomly scattered in the plane. In this case, the caustic also disappears, although 
this happens much more slowly than when Ai -I- A2 > 0. 

A useful quantitative understanding of the steady state can be gained by considering the 
probability distribution of the particle density, P{p). The particle density in the vicinity of a 
caustic (at xq, aligned with the y axis, say) has a divergence of the form p{x, y) ~ {x — xo)~^^^. 
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This indicates that P{p) ^ p^^ for sufficiently large p. However, numerical evaluations of 
P{p) indicates that P{p) ^ p^^ over a significant range of values, only conforming to the 
p~^ prediction for very large p. This surprising observation is a consequence of the nature 
of the folding process which generates the caustics, and we give a brief heuristic argument 
which illustrates the principle. Consider the schematic diagram in Fig. Ol^b), which shows 
a simplified model for the evolution of a fold which forms at f = 0, with a depth which is 
proportional to time t. The corresponding density distribution at time t is P{p,t) = tf{pt) 
for some function /. The pattern contains folds at different stages of their evolution, and the 
observed probability distribution function P{p) is assumed to be a time-average of P{p,t), 
between t = and some cutoff (which can be interpreted as the time between the formation 
of folds). According to this model, 

P{p) - i fdt Pip, i) = ^ r'dx xfix) ^ p-' (6) 
to Jo P to Jo 

(we assume that the integral of xf(x) converges to a constant as its upper limit approaches 
infinity). The model is not applicable at very short times, where the fold resembles the generic 
type illustrated in Fig. 2(a). Fig. |3{c) illustrates the crossover from to p~^ decay in a 
numerical evaluation of Pip)- 

We have argued that the morphology of the steady state caustic patterns can be understood 
qualitatively if we calculate three parameters, namely the rate J at which particles cross 
caustic lines, and the Lyapunov exponents Ai and A2. We will now show how these three 
parameters may be determined quantitatively. The most complete results are available in the 
limit where the correlation time r is very short. 

We linearise eq. (1), and obtain Sr = Sp/m, Sp — —jSp + F(t)(5r, where F is a matrix 
with elements Fijit) = dfiir{t),t)/drj. Consider now the motion of three particles, one 
reference particle and two nearby ones, separated by iSr, Sp) and iSr', Sp') from the reference 
trajectory. The angle Sip separating the vectors Sr and Sr' is very small, and the lengths of 
these vectors, X and X' respectively, are initially equal. We also assume that the vectors 
Sp and Sp' are initially separated by a small angle, of order Sip. We now make a change of 
coordinates from Sr, Sp, Sr', Sp' to the set X, X', 9, Sip, Yi, Y2, Z\, 

Sr = Xns , Sr' = X'he+s^ 
Sp = X{Ying + ¥2-0,0+^/2) , Sp' = X'[iYi + ZiStp)he+8^ + {Y2 + Z2Sip)ng+^+s^] (7) 

ihg denotes a unit vector in two dimensions with direction 9) . We expect that X may increase 
or decrease, Sip decreases with probability one (because random linear mapping of any two 
vectors results in vectors becoming aligned), X'/X remains close to unity, 9 will become 
uniformly distributed on [0, 2tt], and that Yi, Y2, Z\ and Z2 approach a stationary distribution. 
The length of a vector separating two nearby points is X, so that Ai = {X /X). The area 
spanned by the two vectors is 5 A ^ XX' Sip ^ X'^Sip, so that A2 + Ai = (Sip/ip) + 2{X / X) 
(here we used the fact that Sip ^ 1). 

We find the equations of motion for these variables. For X and 5(p we find X — YiX/m 
and Stp = Z2Sip/'m, so that the Lyapunov exponents are 

Ai = {Yi}/m , A2 - Ai + {Z2)/m . (8) 

For 9, we find 9 — Y2/m, and conclude that 9 executes a random walk, becoming uniform 
on [0, 27r] as expected. For the remaining variables we find the following equations, where 
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Yi = -jY, + {Y^-Y,')/m + Fi, 
Zi = -'yZi-2{ZiZ2/2 + YiZi-Y2Z2)/r 



(9) 



F[2 



-7Z2 - 2{Zl/2 + Y1Z2 + Y2Zi)/m - Fi'i 



-P22 



In the limit where the correlation time of the random velocity field is short, we can approximate 
eqns. Q by a set of coupled Langevin equations. We write these in a dimensionless form by 
introducing dimensionless variables Xi, such that (Yi, Y2, Zi, Z2) = m'j{xi, X2, ^3, X4) = mjx, 
t' = jt and find 

dx ^ v{x)dt' + dw (10) 



where vi = —xi + [x"^ 



-X4- 



-X4 - 



a;?), V2 



D, 



'2{xiXi+X2X3), (dwi) 
( 



= ~X2 - 2x1X2-, W3 = -^z - '^■i^A - 2(a;ia;3 - a;2a;4), W4 = 
0, (dwjiduij) = 2Dijdt'. The diffusion matrix with elements 



V -(r+i)/2 





r 

(r+i)/2 






(r+i)/2 
r+i 




'(r+i)/2 \ 


r+i 



where 



Dn 



2m?^^ 



dt (Fn(t)Fn(O)) , T 



1 + 3^2 



(11) 



(12) 



In our simulations, L'o = (3 + a^)(T^/(2TO^7^^^). The three parameters describing caustic 
formation are obtained from the stationary state of the Langevin process (|10f) : the rate J 
at which a representative particle passes through a caustic is the same as the frequency with 
which the area of the parallelogram spanning the two vectors r and r' becomes equal to zero, 
or equivalently the rate at which 5lp passes through zero. An equivalent condition is that 
the trajectory in the {Zi, Z2) plane goes to infinity (reappearing from the reflected direction): 
J = ^j, where j is the rate at which x^ escapes to infinity. Fig. Et^d) shows a comparison 
between a direct evaluation of the Lyapunov exponents (using a method described in [11]), 
and a simulation using eq. H10() . 

Eq. H1U|I is equivalent to a Fokker-Planck equation for the distribution P{x,t), namely 
dP/dt' — V .[—vP + DVP]. Considering a steady-state solution in the limit Dq 0, we 
find that a WKB ansatz, of the form P{x) = exp[~$(a;)/_Do], is appropriate. We therefore 
expect that the escape current vanishes exponentially as Dq 0, being of the form j ~ 
jo exp(— ^oZ-Do) (where jo may have an algebraic dependence on Do). Fig- Ele) shows that 
the caustic formation rate does vanish exponentially as Dq — > 0, with action $0 ~ 0.14. 
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